Goals

The Sulfate Data

The sulfate data in spmodel contains data on 197 sulfate measurements in the continental United States

head(sulfate)

Simple feature collection with 6 features and 1 field
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: 162932.8 ymin: 1080571 xmax: 914593.6 ymax: 1453636
Projected CRS: NAD83 / Conus Albers
  sulfate                 geometry
1  12.925 POINT (817738.8 1080571)
2  20.170 POINT (914593.6 1295545)
3  16.822 POINT (359574.1 1178228)
4  16.227 POINT (265331.9 1239089)
5   7.858 POINT (304528.8 1453636)
6  15.358 POINT (162932.8 1451625)

ggplot(sulfate, aes(color = sulfate)) +
  geom_sf(size = 3) +
  scale_color_viridis_c(limits = c(0, 45)) +
  theme_gray(base_size = 18)

Fitting a Model

We fit and summarize a spatial linear model with an intercept by running

spmod <- splm(sulfate ~ 1, data = sulfate, spcov_type = "exponential")
summary(spmod)

Call:
splm(formula = sulfate ~ 1, data = sulfate, spcov_type = "exponential")

Residuals:
   Min     1Q Median     3Q    Max 
-5.738 -2.605  4.900 13.323 38.099 

Coefficients (fixed):
            Estimate Std. Error z value Pr(>|z|)
(Intercept)    5.924      6.529   0.907    0.364

Coefficients (exponential spatial covariance):
       de        ie     range 
     85.8      10.4 3105165.7 

The broom Functions

Tidy the fixed effect output

tidy(spmod)

# A tibble: 1 × 5
  term        estimate std.error statistic p.value
  <chr>          <dbl>     <dbl>     <dbl>   <dbl>
1 (Intercept)     5.92      6.53     0.907   0.364

Glance at the model fit

glance(spmod)

# A tibble: 1 × 9
      n     p  npar value   AIC  AICc logLik deviance pseudo.r.squared
  <int> <dbl> <int> <dbl> <dbl> <dbl>  <dbl>    <dbl>            <dbl>
1   197     1     3 1140. 1146. 1146.  -570.     196.                0

Augment the data with model diagnostics

augment(spmod)

Prediction (i.e., Kriging)

predict(spmod, newdata = sulfate_preds)

Augment prediction data

aug_preds <- augment(spmod, newdata = sulfate_preds)
print(aug_preds)

Prediction (i.e., Kriging)

Visualize predictions

ggplot(aug_preds, aes(color = .fitted)) +
  geom_sf(size = 3) +
  scale_color_viridis_c(limits = c(0, 45)) +
  theme_gray(base_size = 18)

Your Turn

Your Turn

Goals

Independent Random Error Model

\[ \mathbf{y} = \mathbf{X} \boldsymbol{\beta} + \boldsymbol{\epsilon}, \qquad(1)\]

• \(\boldsymbol{\epsilon}\) is a column vector of random errors.
  • \(\text{E}(\boldsymbol{\epsilon}) = \mathbf{0}\)
  • \(\text{Cov}(\boldsymbol{\epsilon}) = \sigma^2_{ie} \mathbf{I}\)
• \(\mathbf{y}\) is a column vector of response variables, \(\mathbf{X}\) is a design (model) matrix of explanatory variables, and \(\boldsymbol{\beta}\) is a column vector of fixed effects.

Spatial Linear Model

\[ \mathbf{y} = \mathbf{X} \boldsymbol{\beta} + \boldsymbol{\tau} + \boldsymbol{\epsilon}, \qquad(2)\]

• \(\boldsymbol{\tau}\) is a random error vector independent of \(\boldsymbol{\epsilon}\)
  • \(\text{E}(\boldsymbol{\tau}) = \mathbf{0}\)
  • \(\text{Cov}(\boldsymbol{\tau}) = \sigma^2_{de} \mathbf{R}\)
• \(\text{Cov}(\mathbf{y}) \equiv \boldsymbol{\Sigma} = \sigma^2_{de} \mathbf{R} + \sigma^2_{ie} \mathbf{I}\)

Point-Referenced Data

• Data are point-referenced when the elements in \(\mathbf{y}\) are observed at point-locations indexed by x-coordinates and y-coordinates on a spatially continuous surface with an infinite number of locations
  • Consider sampling soil at any point-location in a field
• Distance-based correlation for point-referenced data:
  • Exponential Correlation: \(\mathbf{R} = \exp(-\mathbf{H} / \phi)\)
  • \(\mathbf{H}\) is a matrix of distances, \(\phi\) is the range parameter
• Fit in spmodel using splm() (the workshop’s focus)

Areal Data

• Data are areal when the elements in \(\mathbf{y}\) are observed as part of a finite network of polygons whose connections are indexed by a neighborhood structure
  • The polygons may represent states in a country who are neighbors if they share at least one boundary.
• Neighborhood-based for areal data:
  • SAR Correlation: \(\mathbf{R} = [(\mathbf{I} - \phi \mathbf{W}) (\mathbf{I} - \phi \mathbf{W}^\top)]^{-1}\)
  • \(\mathbf{W}\) is a weight matrix that describes the neighborhood structure, \(\phi\) is the range parameter
• Fit in spmodel using spautor()

Your Turn

Your Turn

Goals

The Moss Data

The moss data contains a variable for log Zinc concentration for moss samples collected near a mining road in Alaska.

ggplot(moss, aes(color = log_Zn)) +
  geom_sf(size = 2) +
  scale_color_viridis_c() +
  scale_x_continuous(breaks = seq(-163, -164, length.out = 2)) +
  theme_gray(base_size = 18)

The Moss Data

The splm() function

The splm() function shares syntactic structure with the lm() function and generally requires at least three arguments

Fit a Spatial Linear Model

Estimation Methods:

• Restricted maximum likelihood (default; likelihood-based)
• Maximum likelihood (likelihood-based)
• Composite likelihood (semivariogram-based) (Curriero and Lele 1999)
• Weighted least squares (semivariogram-based; see weights) (Cressie 1985)

spmod <- splm(formula = log_Zn ~ log_dist2road, data = moss,
              spcov_type = "exponential")
summary(spmod)

The Spatial Linear Model

\[ \mathbf{y} = \mathbf{X} \boldsymbol{\beta} + \boldsymbol{\tau} + \boldsymbol{\epsilon}, \qquad(3)\]

• \(\text{Cov}(\boldsymbol{\tau}) = \sigma^2_{de} \mathbf{R}\)
• \(\text{Cov}(\mathbf{y}) \equiv \boldsymbol{\Sigma} = \sigma^2_{de} \mathbf{R} + \sigma^2_{ie} \mathbf{I}\)
• Exponential Correlation: \(\mathbf{R} = \exp(-\mathbf{H} / \phi)\)

Your Turn

Your Turn

Model Fit

The glance() function returns columns with the sample size (n), number of fixed effects (p), number of estimated covariance parameters (npar), optimization criteria minimum (value), AIC (AIC), AICc (AICc), log-likelihood (loglik), deviance (deviance), and pseudo R-squared (pseudo.r.squared)

glance(spmod)

# A tibble: 1 × 9
      n     p  npar value   AIC  AICc logLik deviance pseudo.r.squared
  <int> <dbl> <int> <dbl> <dbl> <dbl>  <dbl>    <dbl>            <dbl>
1   365     2     3  367.  373.  373.  -184.      363            0.683

Your Turn

Model Diagnostics

The augment() function returns columns with

• .fitted, the fitted value, calculated from the estimated fixed effects in the model
• .hat, the Mahalanobis distance, a metric of leverage
• .cooksd, the Cook’s distance, a metric of influence
• .std.resid, the standardized residual

augment(spmod)

Your Turn

Your Turn

Goals

The Moose Data

The moose data contains moose counts and moose presence for 218 spatial locations in Alaska.

ggplot(data = moose, aes(colour = count)) +
  geom_sf(size = 2) +
  scale_colour_viridis_c(limits = c(0, 40)) +
  theme_minimal(base_size = 18)

The Moose Data

The Moose Prediction Data

The moose_preds data contains spatial locations that were not surveyed.

Fit a Spatial Linear Model

moosemod <- splm(count ~ elev * strat, data = moose,
                  spcov_type = "spherical")
tidy(moosemod)

# A tibble: 4 × 5
  term        estimate std.error statistic p.value
  <chr>          <dbl>     <dbl>     <dbl>   <dbl>
1 (Intercept)   0.310    9.02       0.0344 0.973  
2 elev          0.0141   0.00806    1.76   0.0792 
3 stratM        6.93     2.26       3.07   0.00217
4 elev:stratM  -0.0273   0.0130    -2.10   0.0357 

Predictions

Using predict()

predict(moosemod, newdata = moose_preds)

Using augment()

moose_aug <- augment(moosemod, newdata = moose_preds)
moose_aug

Your Turn

Cross Validation

The loocv() function can be used to perform leave-one-out cross validation on a fitted model object using mean-squared-prediction error (MSPE) loss:

loocv(moosemod)

[1] 32.15933

Your Turn

Goals

Multiple Models

Provide a vector of spcov_types:

spmods <- splm(formula = log_Zn ~ log_dist2road, data = moss,
              spcov_type = c("exponential", "gaussian"))
names(spmods)

[1] "exponential" "gaussian"   

Natural to combine with glances() and predict():

glances(spmods)

# A tibble: 2 × 10
  model         n     p  npar value   AIC  AICc logLik deviance pseudo.r.squared
  <chr>     <int> <dbl> <int> <dbl> <dbl> <dbl>  <dbl>    <dbl>            <dbl>
1 exponent…   365     2     3  367.  373.  373.  -184.     363             0.683
2 gaussian    365     2     3  435.  441.  441.  -218.     363.            0.686

Your Turn

Non-Spatial Random Effects

The random Argument

In splm(), the random argument follows similar syntax to how random effects are specified in the nlme and lme4 packages.

• random = ~ group and random = (1 | group) specify random intercepts for each level of group.
• random = (x | group) specifies random intercepts for group and for each level of group to have a different slope for x.

Non-Spatial Random Effects

randint <- splm(log_Zn ~ log_dist2road,
                data = moss, spcov_type = "exponential",
                random = ~ (1 | sample))

Non-Spatial Random Effects

summary(randint)

Call:
splm(formula = log_Zn ~ log_dist2road, data = moss, spcov_type = "exponential", 
    random = ~(1 | sample))

Residuals:
    Min      1Q  Median      3Q     Max 
-2.6234 -1.3228 -0.8026 -0.2642  1.0998 

Coefficients (fixed):
              Estimate Std. Error z value Pr(>|z|)    
(Intercept)    9.66066    0.26770   36.09   <2e-16 ***
log_dist2road -0.55028    0.02071  -26.58   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Pseudo R-squared: 0.6605

Coefficients (exponential spatial covariance):
       de        ie     range 
3.153e-01 2.094e-02 1.083e+04 

Coefficients (random effects):
1 | sample 
   0.07995 

More Complex Random Effects Models

yearint <- splm(log_Zn ~ log_dist2road,
                      data = moss, spcov_type = "exponential",
                      random = ~ (1 | sample + year))
yearsl <- splm(log_Zn ~ log_dist2road,
                      data = moss, spcov_type = "exponential",
                      random = ~ (1 | sample) + 
                       (log_dist2road | year))

Your Turn

Anisotropy

An anisotropic covariance does not behave similarly in all directions

• Consider a spatial covariance influenced by the prevailing wind direction.

aniso <- splm(log_Zn ~ log_dist2road,
              data = moss, spcov_type = "exponential",
              anisotropy = TRUE)
glances(spmod, aniso)

# A tibble: 2 × 10
  model     n     p  npar value   AIC  AICc logLik deviance pseudo.r.squared
  <chr> <int> <dbl> <int> <dbl> <dbl> <dbl>  <dbl>    <dbl>            <dbl>
1 aniso   365     2     5  362.  372.  372.  -181.     363.            0.705
2 spmod   365     2     3  367.  373.  373.  -184.     363             0.683

Anisotropy

summary(aniso)

Call:
splm(formula = log_Zn ~ log_dist2road, data = moss, spcov_type = "exponential", 
    anisotropy = TRUE)

Residuals:
    Min      1Q  Median      3Q     Max 
-2.5279 -1.2239 -0.7202 -0.1921  1.1659 

Coefficients (fixed):
              Estimate Std. Error z value Pr(>|z|)    
(Intercept)    9.54798    0.22291   42.83   <2e-16 ***
log_dist2road -0.54601    0.01855  -29.44   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Pseudo R-squared: 0.7048

Coefficients (exponential spatial covariance):
       de        ie     range    rotate     scale 
3.561e-01 6.812e-02 8.732e+03 2.435e+00 4.753e-01 
attr(,"class")
[1] "exponential"

Your Turn

Your Turn

Partition Factors

A factor variable (i.e., categorical) is a partition factor when two observations in separate levels of the partition factor should be uncorrelated

• Observations from the same year are spatially correlated but observations from different years are not:

part <- splm(log_Zn ~ log_dist2road,
             data = moss, spcov_type = "exponential",
             partition_factor = ~ year)

Fixing Covariance Parameters

Steps:

1. Use spcov_initial() to specify the covariance type and any known, fixed covariance parameters.

2. Instead of specifying the spcov_type argument in splm(), specify the spcov_initial argument with the object from (1).

init_spher <- spcov_initial("spherical", range = 20000, known = "range")
splm(log_Zn ~ log_dist2road, data = moss,
     spcov_initial = init_spher)

Your Turn

Your Turn

Random Forest Spatial Residual Models

See workbook for an example of prediction using random forest spatial residual models

Goals

Large Data Challenges

• Inversion of \(\boldsymbol{\Sigma}\) for data sets with around 10,000 or more observations is challenging (and often infeasible) on a standard computer
• spmodel implements “local” spatial indexing (Ver Hoef, Dumelle, et al. 2023) for model fitting
  • Induce some sparsity in the covariance matrix
  • Set local = TRUE in splm()
• See workbook for example
  • 5,000 observations fit in 12 seconds (no parallel)

Large Data Challenges

• spmodel uses “local neighborhood prediction” to predict for unobserved spatial locations for a model fit to a large data set
  • Only a subset of observations are used to predict the response at a particular location
  • Set local = TRUE in predict() or augment()
• See workbook for example
  • 3,000 predictions in 11.7 seconds (no parallel)

Goals

The Spatial Generalized Linear Model

\[ g(\boldsymbol{\mu}) = \mathbf{X} \boldsymbol{\beta} + \boldsymbol{\tau} + \boldsymbol{\epsilon}, \]

• \(g(\boldsymbol{\mu})\) is the link function that “links” a function of the mean of \(\mathbf{y}\) to \(\mathbf{X} \boldsymbol{\beta}\), \(\boldsymbol{\tau}\), and \(\boldsymbol{\epsilon}\)
• Fit models using a novel application of the the Laplace approximation (Ver Hoef, Blagg, et al. 2023)

Response Distributions

The Moose Data

Model Fitting

poismod <- spglm(count ~ elev * strat, data = moose,
                 family = poisson, spcov_type = "spherical")

• The family argument can be binomial, beta, poisson, nbinomial, Gamma, or inverse.gaussian

• Notice the similarities between spglm() and glm()

Model Fitting

summary(poismod)

Call:
spglm(formula = count ~ elev * strat, family = poisson, data = moose, 
    spcov_type = "spherical")

Deviance Residuals:
    Min      1Q  Median      3Q     Max 
-1.4245 -0.7783 -0.3653  0.1531  0.5900 

Coefficients (fixed):
             Estimate Std. Error z value Pr(>|z|)    
(Intercept) -2.230575   0.958201  -2.328 0.019919 *  
elev         0.007623   0.003129   2.437 0.014820 *  
stratM       2.752234   0.782853   3.516 0.000439 ***
elev:stratM -0.010248   0.004472  -2.292 0.021928 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Pseudo R-squared: 0.09573

Coefficients (spherical spatial covariance):
       de        ie     range 
    3.892     1.163 51204.657 

Coefficients (Dispersion for poisson family):
dispersion 
         1 

Your Turn

Prediction

Using predict()

predict(poismod, newdata = moose_preds)

Using augment()

augment(poismod, newdata = moose_preds)

Additional Modeling Features

All advanced features available in spmodel for spatial linear models (splm()) are also available for spatial generalized linear models (spglm())

• Fit and predict for multiple models
• Non-spatial random effects
• Anisotropy
• Etc.

Your Turn

Goals

Simulating Gaussian Data

1. Use spcov_params() to specify the correlation structure and covariance parameters

params <- spcov_params("exponential", de = 1, ie = 0.5, range = 5e5)

2. Use sprnorm(params, data) to simulate a realization of the spatial process defined by spcov_params() at locations in data.

set.seed(1)
sulfate$z <- sprnorm(params, data = sulfate)

Simulating Gaussian Data

3. Visualize

ggplot(sulfate, aes(color = z)) +
  geom_sf(size = 2) +
  scale_color_viridis_c() +
  theme_gray(base_size = 14)

The Empirical Semivariogram

What does an empirical semivariogram of the simulated data look like?

esv_out <- esv(z ~ 1, sulfate)
ggplot(esv_out, aes(x = dist, y = gamma, size = np)) +
  geom_point() +
  lims(y = c(0, NA)) +
  theme_gray(base_size = 14)

Simulating Other Data

Use sprpois(params, data) to simulate a Poisson realization and visualize

sulfate$p <- sprpois(params, data = sulfate)
ggplot(sulfate, aes(color = p)) +
  geom_sf(size = 2) +
  scale_color_viridis_c() +
  theme_gray(base_size = 14)

A Caution

• Simulating spatial data in spmodel requires the Cholesky decomposition of the covariance matrix, which can take awhile for sample sizes exceeding 10,000
• Regardless of the number of realizations simulated, this Cholesky decompsition is only needed once
• This means that simulating many realizations (via samples) takes nearly the same time as simulating just one.

Goals

The Seal Data

ggplot(seal, aes(fill = log_trend)) +
  geom_sf() +
  scale_fill_viridis_c() +
  theme_bw(base_size = 18) 

Figure 8: Distribution of the seal data.

Model Syntax

Model syntax for spautor() (spgautor()) is similar to the syntax used for splm() (spglm()):

sealmod <- spautor(log_trend ~ 1, data = seal, spcov_type = "car")

Model Output Interpretation

summary(sealmod)

Call:
spautor(formula = log_trend ~ 1, data = seal, spcov_type = "car")

Residuals:
     Min       1Q   Median       3Q      Max 
-0.34441 -0.10403  0.04423  0.07351  0.20489 

Coefficients (fixed):
            Estimate Std. Error z value Pr(>|z|)   
(Intercept) -0.07103    0.02492  -2.851  0.00436 **
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Coefficients (car spatial covariance):
     de   range   extra 
0.03226 0.42020 0.02235 

Your Turn

Your Turn

Prediction

Predictions must occur for “missing” (NA) responses in the observed data

Using predict()

predict(sealmod)

Using augment()

augment(sealmod, newdata = sealmod$newdata)

Your Turn

is_missing <- is.na(seal$log_trend)
seal_nomiss <- seal[!is_missing, , ]

Thank You!

• Thank you so much for attending
• Please reach out with comments / questions / suggestions / bugs (Dumelle.Michael@epa.gov)
• A6.02 16:00 - 16:15 Thursday, 20 July, 2023 UMC Glenn Miller Ballroom Middle
  • Slides available for download https://github.com/USEPA/spmodel.spatialstat2023/blob/main/spmodel-slides-intro.html